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APPLICATION OF LARGE-EDDY INTERACTION MODEL TO CHANNEL FLOW 
S. K. Hong* and M. W. Rubesin** 

NASA Ames Research Center 
Moffett Field, CA 94035 

A procedure utilizing an expansion of proper orthogonal functions (or modes) 
to predict a fully developed flow in channel is derived. To examine numerical 
and conceptual difficulties, preliminary computations are performed with 
assigned mean velocity, and turbulence is expressed with only the first 
mode. The nonlinear interactions of the components of the first mode are 
treated specifically, with the influence of higher modes neglected; this 
treatment requires adjustment of the skewness and effective Reynolds number to 
assure energy equilibrium of the first mode. Computational results show that 
the first mode possesses the structural character similar to that of the 
entire flow. 


Nomenclature 

d = channel half-width 

H = channel width 

k = ( k -| , k ^ ) > wave number vector 

* -fi * k 3 

p = pressure fluctuation 

Pi , . . . , P q = large-eddy spectra 


2 

q 

Re 


R ij 

s 

t 

u. 


2 2 2 

= u + v + w , twice of turbulent kinetic energy 

= Reynolds number 

= two-point velocity correlation 
3 2 3/2 

= , skewness parameter 

= tU Q /H, dimensionless time 
= instantaneous velocity 
= time-averaged velocity 
= mean velocity at the channel centerline 
= (u,v,w), velocity fluctuation 
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= wall-friction velocity 
= Reynolds stress tensor 


u* 

u. u . 

_i J 

2 

u = streamwise turbulent intensity 

uv = shear stress 

x = (x,y,z), position vector 

Greek Symbols 

a = random coefficient 

n _ 

A =Va 1 2 | energy of nth mode 

n r n 

v = kinematic viscosity 

v, p = eddy viscosity 

= nth mode of pressure fluctuation 
= nth mode of velocity fluctuation 

Subscripts 

i = 1,2, or 3 in streamwise, normal, or spanwise direction, 

respectively 

n = indicates order of mode 

Superscripts 

1 = indicates first mode 

n = indicates nth mode 

0 = ensemble average 

( )' = dif ferentiation with respect to y 

1 Introduction 

The mean velocities and the local structure of turbulence under given 
initial and boundary conditions are subjects of importance in the engineering 
predictions of inhomogeneous, turbulent shear flows. The most popular predic- 
tion methods assume that the Navier-Stokes equations are adequate for describ- 
ing turbulent flow on an instantaneous basis and are used to develop statisti- 
cal equations for the various turbulent moments , including the Reynolds 
stresses. However, the Navier-Stokes equations involve more turbulent moments 
than there are equations to explain them, forming an unclosed system. Various 
levels of closure schemes have been proposed: the zero-equation model [1], 

the two-equation model [2,31, and the Reynolds stress equation model C^,5]. 
Each of the foregoing methods requires the introduction of several empirical 
constants with respect to various turbulent processes, and only provides 
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approximate predictions of the nature of individual turbulent processes aris- 
ing in a given flow. An approach that examines the dynamics of the turbulence 
may require less reliance on modeling. Earlier, one of the authors used the 
Proper Orthogonal Decomposition Theorem (PODT) [6,7] to analyze the structural 
character of turbulence in boundary layers over variously curved walls [8]. 

He found that the first mode of the decomposition exhibits the structural 
character of the averaged turbulence moments quite well even though all the 
nonlinear effects were modeled in the framework called the Large-Eddy Interac- 
tion Model (LEIM) [9~11]. This work suggested that only a few properly iden- 
tified modes in the decomposition may be needed if the PODT is utilized in a 
predictive sense and led to the current work. 

It should be noted that Lumley’s [7] primary purpose was to define unam- 
biguously the meaning of an "eddy." Given two-point velocity correlations 
from experiments, the PODT has been applied to a pipe flow [12], a wake [13], 
and a flat-plate boundary layer [14] as a means of extracting the features of 
a dominant eddy. In a similar approach with the results from a computation- 
ally simulated, fully developed channel flow, Moin [15] has specifically 
investigated as to how many modes are necessary in the PODT to reproduce the 
turbulent intensities and shear stress. Moin [15] shows that in the case of 
shear stress it takes the first 15 terms in the series representation before 
the calculated stress distribution matches the value simulated earlier [16] 
across the boundary layer. However, when only the wall region is examined, 
the sum of the first three modes yields the "experimental" results quite 
well. This small number of modes needed to describe the turbulence gave 
further encouragement to the authors to extend the PODT approach as a predic- 
tive tool as was done in the LEIM. 

The basic procedure of the LEIM consists of the following steps: 

(a) decomposing the velocity fluctuations into orthogonal functions with 
random coefficients; (b) constructing dynamical equations for those functions; 
(c) identifying the first mode as an organized structure that contributes most 
to the energy [7,17]; and (d) evaluating the large eddy which interacts with 
the mean flow and the eddy-eddy interactions. 

In the past, all the nonlinear terms in the LEIM were modeled in a linear 
form utilizing either an anisotropic eddy viscosity or a diffusion velocity. 

In this process three empirical constants were introduced in the closure and 
then determined by matching shapes between normalized Reynolds stresses 
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calculated from the first mode and experimental measurements. In view of the 
emphasis of the past applications of the LEIM on evaluating the normalized 
structure of the first mode, it was primarily a diagnostic method. In the 
present work to develop a predictive method, the turbulent transport processes 
have been reexamined and retained in their nonlinear form to minimize the 
dependence on turbulence modeling and to allow evaluation of the magnitudes of 
the moments. The applicability of the new transport model has been illus- 
trated in a channel flow that is inhomogeneous in the direction normal to the 
wall. The computed results shown here are restricted to use of only a single 
mode in the decomposition and, as stated earlier, this forced the introduction 
of the skewness and the effective Reynolds number as parameters of the prob- 
lem. Although the current work has not yet demonstrated the uniqueness of a 
set of these parameters, results based on different sets of these parameters 
that satisfy the energy equilibrium of the large eddy yield Reynolds stresses 
that are less than a few percent apart over the entire channel. 


2 Large-Eddy Interaction Model 

One can consider a decomposition of the instantaneous velocity, l^, into 
a time-mean value and it3 fluctuation as 


U. = U + u. 
i i 1 


( 1 ) 


where is the time-averaged mean velocity. The velocity fluctuation in 

equation (1) may be represented as a series in terms of orthogonal functions, 
n = 1,2,3-..}, that is 


u i (x,t) = £ a (t)i/ n ^(x) 
n=1 


(2) 


Here a (t) are random coefficients and uncorrelated from one another, which 
n 

means 


a a = A (t ) • 6 
m n n mn 


(3) 


The overbar, ( ) > represents an ensemble average. In addition, the 


are assumed to be orthonormal functions: 


S' 


, (p Y (q) dx = 6 

i pq 


(4) 
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It then follows [7] that the two-point velocity correlation, R n , can be 
expressed with X^ and <j>^ in the form: 


R i j(x,x",t) » u i (x,t)Uj (x",t) 


- E X^{ n) (x) • <f> (n) (x") (5) 

n=1 J 


Equation (5) implies that if one can predict (A then various two-point 

correlations and Reynolds stress components can be calculated directly. 

To formulate a framework for we first substitute equation (2) 

into the Navier-Stokes equation for u^ . After a few manipulations, dynamical 
equations for the nth mode then become the following under incompressible 
flow assumptions [7]: 


9U n <! > i (n) ) 

at 


* u j ^ <vi n>) 


3U. , , 

(X <j> n) ) 

3Xj n Y j 


_3 

3x 


{ » “ a a a . s ^ 

E E -V- 2 * p V, q) [ 

p = 1 q=1 n J J 


3( X Tt (n) ) 
n 

3x, 


+ v 


3 2 (X n <0i (n) ) 

3x. 3x. 

J J 


(n 


The continuity relationship yields 




3x. 


- = 0 


1,2,3, ...) 


(sum over index j ) 


( 6 ) 

(7) 


where (X^ir^) = -a n p/p. However, only the system of equations corresponding 
to n = 1 in the above will be considered in the current paper. The index, 
i, may have values of 1, 2, or 3, corresponding to the streamwise direction 
(x), the local normal to the wall (y), and the spanwise direction (z), respec- 
tively. The nonlinear term appearing in equation (6) when n = 1 represents 

the interaction process between eddies of various order and requires a closure 

(2) (3) 

assumption in view of the presence of higher order modes, <j>. ^ ,.... 
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2.1 Transport Process 

The nonlinear interaction between the first mode is retained in its 
original form, but the interactions involving modes higher than are 

modeled. A simple way of accounting for the effect of the higher modes is to 
group them together and to relate this effect to a known quantity. An eddy 
viscosity is introduced for this purpose [18,19]: 


00 00 


£ E 

p = 1 q=1 


a, ct ct / \ / \ 
i J 



(1) 
i ' 


( 1 ) 


X 3x j 


) 


+ 


3 ( X , j 


( 1 ) 


9x. 

i 




( 8 ) 


which has also been suggested by Lumley [7]. In equation (8), denotes an 
eddy viscosity. In the present analysis, v T is kept equal to a constant, 
independent of mesh dimensions or the distance from the wall. This has an 
effect of reducing the effective Reynolds number by a factor of v/(v + v^). 
The purpose of the current study is to investigate whether we need to include 
the effects of higher modes. Although improved models (with variation in 
the y-direction) are required, in the present preliminary study a constant 
eddy-viscosity model is incorporated to achieve a steady-state solution. We 
realize that the computational results indicating the magnitude of the effect 
of higher modes will be affected by this choice. 

Upon substituting equation (8) into (6), one obtains a closed system of 
equations for 


J J 


( 1 ). 


3U. 


afr <v (1) > *<»♦ Vjrsr'Vi "’ 1 

i J J 


(9) 


3 2 3/2 

where S (= 0 ^/( 0 ^) ) is the skewness factor of the random coefficient, 

a 1 . Assuming that the mean velocity is given, the system of equations 
involves then a structural parameter, S, and a stability parameter, v^, which 
need to be chosen. For the LEIM to be a completely predictive scheme, it 
requires incorporation of the mean momentum equation for along with 
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equations (7) and (9). However, for the present, emphasis will be placed on 
how the large eddy interacts and reacts to a known mean flow field. 

3 Application 

Fully developed turbulent channel flow is a case for which a large data 
base has become available over the years. Accordingly, two-dimensional 
channel flow has been chosen to demonstrate the validity of the approximation 
for the nonlinear eddy-eddy interaction terms as proposed in equation (8) and 
to establish the contribution of the first mode to various statistical turbu- 
lence quantities. In the fully developed region of the channel, the mean 
velocity is one-dimensional and is dependent only on the normal coordinate, y, 
where y = 0 corresponds to the lower wall and y = H to the upper wall. 
Thus, the turbulent flow in the two-dimensional channel flow can be regarded 
as homogeneous in both streamwise (x) and spanwise (z) coordinates, while 
strong inhomogeneity is retained in the y-direction. For this case, one can 

define spectral functions for the first mode of velocity fluctuation, X.<J>f^, 

( 1 ) 1 

and of pressure fluctuation, X^tt , as follows. 

<t>,0< 1 ,y,k ,t) = — — - IT {X (t)<t>, (1) (x,y,z)} 

1 1 3 (2ir) 2 *'•'-» 1 1 

x exp{-i(k.jX + k^z)}dx dz 

( 1 °) 

ir(k ,y,k t) = — ~ if {X. (t)ir^ 1 ^(x,y,z) } 

1 3 (2tt) 2 1 

A 

x expf-iCk.jX + k^z)}dx dz 

a 

where k 1 and k^ are wave numbers and i = /-I . The superscript indicating 
mode (1) is understood in the spectral functions. Applying these definitions 
into equations (7) and (9), one obtains four complex equations with respect to 

A A 

the large eddy spectra, (i = 1,2,3) and m. These equations can be further 
divided into eight equations for the real and imaginary parts defined accord- 
ing to the following notation. 
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( 11 ) 


♦l = P 1 + 1P 2 


*2 “ P 3 + iP 4 
♦ 3 " P 5 + * P 6 

/S. A 

TT = + iPg 

The nonlinear terms require the convolution theorem [20] during the 
transformation of the system of equations, equations (7) and (9), into the 
mixed, (l<i ,y ,kg,t ) , space. The spectral equations become 

A 

9tr + ik i u<t) i + u ’^2 + s 


34>o 

W * ‘ k 1 U *2 * S • 


- ik^w - (v + v T > C4>” - k 2 <f> 1 ) = 0 


-IT' - (v + v T ) (<i>” - k <J> 2 ) = 0 


9<j>, 

w 


‘ k , u * 3 * s ' ‘♦"’♦j"’} 


A A 


- ikgir - (v + v T )(<(>" - k <|) ) = 0 


( 12 ) 


(13) 


(14) 


ik 1 4> 1 + + ik^ = 0 


(15) 


2 2 2 

where ( )' = d( )/dy and k = k^ + k^. From equations (12)~(14), we denote 



i k * * 4> ^ (k M ) <J> , (k 


k") 


d<J) . (k") 

+ 

dy 


(k - k") 


+ ik”<J. i (k")<(. 3 (k - k") dky dk^ 


( 16 ) 


where k = (k 1 .k^) . 

From a consideration of two-point velocity correlations and Reynolds 
stresses, a relation can be found between the spectra of the double velocity 

A 

correlations and the <Jk's [21]: 
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(17) 


R[] ) (k! | ',y,k”,t)5(k + k") 


<|) i (k 1 ,y,k 3 ,t)<|»j(k^,y,k^,t: 


where 5 is the Dirac delta function and superscript (1) indicates the con- 

a 

tribution of the dominant mode. R. . is defined in 


R iJ (r 1 ,y,r 3 ,t) 


JJ R i j(k 1 ,y,k 3 ,t)exp{i(k 1 r 1 + k^) }dk 1 dk 3 (18) 


For a flow which is homogeneous in (x ,y)-planes : 


R i j(x,y,z,t;r 1 ,r 3 ) 

B u^x.y.z.tjujta + r 1 t y,z + r~£T 

co - 

= 2 0 t^(t)<|>. (n) (x,y ,z)«j> ( . n) (x + r ,y,z + r,) 
n=1 J J 

= o^[ 1) (x,y,z)4)j 1) (x + r 1 ,y ,z + r 3 ) (19) 


where = indicates simple truncation after the first mode. When 

r 1 = r^ = 0 in equations (18) and (19), the two-point correlation reduces to 

the usual Reynolds stress tensor, u^Uj(y,t). Thus, at time t: 

Vj ( 1 ) (y,t) 

A A 

* 1 (k 1 ,y,k 3 ,t)<j>*(k 1 ,y,k 3 ,t)dk 1 dk 3 (20) 

where ( )* denotes the complex conjugate of ( ). The following structural 

2 2 

quantities are also presented: (i) normal stress intensities (u./q ), where 

2 2 2 2 2 1 

q - u + v ♦ w ; (ii) shear stress intensity (uv/q ); (iii) orientation of 

the principal axes of the large eddies (0), which is given by 
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0 - (1/2)tan" 1 {-2uv/(u 2 - v 2 )} 


2 2 

and (iv) anisotropy (u /v ). 


4 Calculation Procedure 

The next objective is to calculate the first mode in the velocity fluctu- 
ation expansion and to investigate the extent of the influence of the first 
mode on the turbulent intensities and the other various structural quanti- 
ties. Attention is given to the role of the skewness parameter, S, that 
appears in equation (9), on the solution for the first mode, ^ , and on the 
structure of turbulence deduced therefrom. 

To solve equations (12)-(15), after one adopts a numerical algorithm he 
or she needs: (i) initial and boundary conditions, (ii) the local mean veloc- 

ity profile, and (iii) a proper skewness factor as well as an eddy viscos- 
ity. Numerical results of the statistical quantities can then be compared 
with the classical measurements owing to Laufer [22], for example, from which 
a particular flow condition is selected as 

U Q * 7.574 (m/sec) 
u# - 0.2891 (m/sec) 

H = 12.7 (cm) 

Re = U H/v - 60380 
o 

where U Q , u*, H, and Re are the mean velocity at the channel centerline, 
wall-friction velocity, channel width, and Reynolds number, respectively. 

4.1 Initial and Boundary Conditions 

A numerical solution for the large eddy spectra governed by the system of 
equations (12)— (15) is determined as an initial-boundary value problem in the 
(y,t) space for various values of wave number, k. The initialization can be 
arbitrary because of the goal of achieving a steady solution in the presence 
of a fixed-mean strain. 

At the wall (y =0), the no-slip condition requires 

P 1 “ P 2 " P 3 " P 4 * P 5 “ P 6 “ 0 (22) 

The boundary condition for the pressure spectrum is deduoed again from the 

v component (normal to surface) equation applied at the wall [23]. The 

spatial derivative for the pressure fluctuation is then discretized using a 

three-point forward formula beginning at y « 0. 
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For the other set of boundary conditions, the flow field has been assumed 
to be symmetric with respect to the centerline (y = H/2) , giving 

p j = P' 2 = = P£ = P} = P£ = 0 (23) 

P 3 = P 4 - 0 (24) 

where ( )' denotes the derivative with respect to y. It should be pointed 
out that the entire channel from y = 0 to y = H has been solved in a 
single instance with the no-slip condition imposed at both ends, y = 0 and 
y = H. The results show symmetric profiles for u and w spectra and antisym- 
metric for v spectra (P^ and P^), justifying the use of the current boundary 
conditions at the centerline. 

The implicit numerical scheme employed [24] utilizes a two-point backward 
differencing in time and a three-point central differencing in y. The non- 
linear convolution integrals are treated as known by evaluating them at a 
previous time step when the solution is known. The numerical integration for 
these terms is carried out employing the trapezoidal rule over the wave number 
space, (k^.kjp, at a given point, (k^y.k^.t). This enables the formation of 
a system of matrix equations for all y at each advanced time level, where 
the coefficient matrix becomes block tridiagonal and diagonally dominant. The 
inhomogeneous y-coordinate is discretized as suggested by Murphy and Rubesin 
[25] and the half of the channel is divided by 35 nonuniform grids. 

The wave number space, (k-j.k^), has also been divided into strongly 
nonuniform meshes. The wave number plane is covered with (17,17) grids where 
the values are equally spaced in the logarithmic scale for each wave number 
direction in the range between -10 and 10 (1/cm). It is found (via numerical 
experimentation) that the wave numbers outside this range have negligible 
effect on the turbulent stresses. 

4.2 Mean Velocity Profile 

The mean velocity profile is approximated by a near-wall Prandtl-Taylor 
model and a blending profile near the center plane. 

(i) U + = y + (y + < 12) 

(ii) U + = 3.0 In y + + 5.5 (12 _< y + < 760) 

(iii) U/U Q = 1.0 + 0.068 log (y/d) (y + >_ 760) 

where U + = U/u*, y + = yu*/v, and d = H/2. The mean velocity profile in the 
outer layer (iii) has been modified in this form for the purpose of matching 
with the law of the wall (ii) smoothly using the experimental data. 
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4.3 Effect of Parameters 

During the computation, a fairly small time step (normalized by the mean 
velocity at the channel centerline and by the channel height) of about 0.001 
has to be used to ensure numerical stability. Although , governed by 
equations (12)— (15), can vary in time in any fashion with the current eddy- 
viscosity model, the fully developed steady-mean flow and boundary conditions 
used here cause the solution to converge to a steady state. In the current 
computation, the procedure has been continued up to 200 iterations in time to 
achieve an accurate asymptotic solution. 

5 Results and Discussion 

The effect of various values of the parameters, S and v^,, in equa- 
tions (12)-(14) on the solution is discussed. After a pair of those values is 
selected, numerical results obtained in the mixed space (k^.y.k^.t) are inte- 
grated over the wave number space to yield Reynolds stresses as a function 
of y. The computed stresses and structural quantities are compared with the 
experimental data of Laufer [22]. 

5.1 Determination of Parameters (S and v T ) 

The skewness parameter, S, can be regarded as a structural parameter 

/\ 

affecting the anisotropy obtained from the solution, <j>^. On the other hand, 
the primary role of the eddy viscosity, v T , is to stabilize the growth of the 
solution subject to production in a fixed-mean velocity field. We pick a 
value of S first and then determine corresponding value of v T which yields 
a steady-state solution. One may argue that a choice of a particular set of 
S and v T is not unique on the ground that other combinations of S and v T 
could also produce steady-state results. It has been found, however, that the 
solution is rather insensitive to the choices of the combinations of S and 
that yield steady solutions. 

Figure 1 shows the growth of the turbulent kinetic energy, integrated 
over the channel, for the various values of S as a function of time in the 
absence of higher modes, v T = 0, and for a single case with v T = 18. With 
v T = 0, the growth rate increases with increasing S, and it was found for 
| S | > 0.03 the solution grows so rapidly that it becomes unstable. However, 
for S = 0.01 and v^, = 18, the desired steady state in kinetic energy is 
achieved. Also, for S = 0.03 a steady solution for kinetic energy occurs 
when v T = 22. 
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To examine further whether each component of the kinetic energy has 
indeed reached a steady state for the above two sets of parameters, each 
component of the three turbulent intensities is integrated over the channel 
from y = 0 to y = d. The results are presented in Fig. 2 as a function of 
time. It turns out that the u-component energy for S = 0.01 and v T = 18 
maintains a constant value, but that for S = 0.03 and v T = 22 continues to 
decrease slightly in Fig. 2. For both cases the w component continues to 
increase whereas the v component decreases, again at a slow rate. The 
behavior of these different growth rates can be attributed to the use of an 
isotropic eddy viscosity in equation (8) and suggests the use of an aniso- 
tropic eddy viscosity or some other alternative. Nevertheless, in view of the 
small growth rate in the v and w components, this behavior is believed to be 
relatively unimportant so that no attempt was made to eliminate this continu- 
ally varying anisotropy. The authors favor the case of S = 0.01 and 
v T = 18 because it yielded a steady solution slightly better than the other 
case. However, the calculated Reynolds stresses for these two sets of param- 
eters were well within a few percent of each other (see Table 1 for numerical 
values of the Reynolds stresses). Thus, the choice of which set to use is not 
critical to the results shown in the rest of the paper which is based on 
S = 0.01 and v T = 18. 

5.2 Eddy-Eddy Interactions 

Let's write the eddy-eddy interactions affecting the net production of 
the first mode in the form: 


rHJ 

'Mf’’ * s if- 

> v. J J 



"B" 


9 2 (A 1 <))[ 1) ) 

y T 9x.3x. 

V 

"C" 


(25) 


The term "A" governs the interaction with the mean flow. The term "B" corre- 
sponds to the transport of the whereas the term "C" represents the 

effects of the higher modes modeled by the eddy viscosity. Figures 3~5 illus- 
trate the nature of interactions among the eddies which are identified with 
"B" and "C" in equation (25). The large-eddy/large-eddy interaction in the 
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Pi equation [which is the real part of equation (12)] is used as an illustra- 
tive example and is shown for three sets of wave numbers, 0.2, 0.5, and 1.0 
(1/cm), in Fig. 3. It is found that the values of "B" in equation (25) for 
wave numbers the same as or less than 0.2 (1/cm) are predominantly negative in 
the inner part of the boundary layer. A negative value of "B" refers to an 
energy supply and a positive value, an energy drain because the value of 
itself in the inner layer is negative for those wave numbers. The profiles of 
the nonlinear term show both types of behavior at the wave number about 0.5 
(1/cm) and positive behavior at higher wave numbers than 0.5 (1/cm). Thus, 
the nonlinear eddy-eddy interactions for lower wave numbers cause energy gain 
while that for higher wave numbers dissipates the energy. 

Shown in Figs. 4 and 5 are comparisons between "B" and "C" terms in 
equation (25) for the P^ equation. Figure 4 shows terms "B" and "C" at the 
low wave number k 1 =0.1 (1/cm) and Fig. 5 compares the same terms at the 
value of k.| =5.0 (1/cm), both for a fixed value of the wave number k^ at 
0.05 (1/cm). Profiles of the P^ are also provided in the two figures to 
indicate their behavior in y at the same wave numbers. An opposite sign in 
"B" or "C" from that of P^ implies energy loss and the same sign implies 
energy gain. The value of k^ = 0.05 (1/cm) was chosen in these illustrations 
because it emphasizes the differences occurring in the alternatives of k ^ . 
Similar results are expected of "B" and "C" for other values of k^ when 
k^ is varied in the same manner. Figure 4 shows that at the same wave number 
(k^ =0.1) and y/d = 0.5, the turbulent energy transfer because of the first 
mode self-interactions is much smaller than that caused by the rest of the 
modes when the latter are modeled with = constant. If other models were 
employed for the eddy-eddy interactions, this emphasis on energy loss could be 
reduced. The magnitude of the higher mode-interactions exceeds that of the 
first mode self-interactions considerably in both Figs. 4 and 5 to dissipate 
the energy gained not only from the large eddy but also from the mean flow 
[term "A," equation (25)]. Again from the signs of P^ , "B", and "C" in 
Figs. 4 and 5, one can infer that the higher modes interactions drain the 
energy consistently for all wave numbers, while large eddy self-interactions 
either supply energy for lower wave numbers or subtract it for higher wave 
numbers. Study of these figures illustrates that the influence of the eddy- 
viscosity model for the higher modes interactions is dominant and suggests 
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multiple mode analysis (in conjunction with improved models for higher modes 
interactions) is necessary to capture the nonlinear process properly. 

5.3 Reynolds Stresses 

From the deterministic large-eddy spectra, through Pg, Reynolds stress 

components have been obtained from equation (20). In Figs. 6 and 7, the 

2 2 

Reynolds stresses obtained as a function of (y,t) are given for u and v 
only at every 25 time steps to show the development in time. The solution 
adjusts itself quickly in time and the effect of the initial conditions 
appears to be minimal. To see whether the Reynolds stresses in Figs. 6 and 7 
change their profiles, the u and v profiles at t = 0.1 5, 0. 175, and 0.2 are 
compared in Fig. 8 as a function of y/d. It shows that the u and v compo- 
nent intensities have indeed achieved an equilibrium profile for t > 0.15. 

The same observation has been made for other Reynolds stress components. 

For detailed comparison, the Reynolds stresses are given as a function 
of y/d at time t = 0.2 in Fig. 9 along with experimental distributions 
taken from Laufer [22]. The first mode contributes approximately 3G£ of the 
observed intensities, although the shape agrees in general trend with the 
experimental distribution. In case of the w component, the use of the 
isotropic viscosity has caused spuriously higher proportion of the calcu- 
lated w component than the u or v component in the contribution to the 
energy . 

5.4 Structural Quantities 

The single-point structural quantities defined in Section 3 are calcu- 
lated from the Reynolds stresses and are compared with corresponding quanti- 
ties obtained from measurements of Laufer [22]. In spite of low intensity 
levels of the large eddy as shown in Fig. 9, the normalized structural quanti- 
ties in Fig. 10(a-d) agree well with the corresponding experimental quantities 
except for the anisotropy of the flow in the outer part of the boundary layer 
shown in Fig. 10(d). It is believed that this latter result is the conse- 
quence of the imposed symmetry at the channel centerline, that is, v = 0. 

6 Concluding Remarks 

Since the current paper deals only with the first mode and the balance 
regarding the first mode's origin, maintenance, and destruction, the question 
as to the number of modes in the series required to represent the turbulence 
field accurately to predict engineering quantities still remains unanswered . 


15 



However, the first mode is shown to be so significant that it supplies about 
30? of turbulent kinetic energy and possesses a structural character that 
closely matches the experimental trends of the entire turbulence field. 

When the turbulent transport is truncated to only the interactions 
between the components of the first mode, it is found that the drain of energy 
through the first mode is insufficient to balance the energy gain from the 
mean motion. Thus, the nonlinear interactions of higher modes are necessary 
in the dynamics of the first mode. However, the constant eddy viscosity used 
in the current closure is inadequate. For the improvement of the closure with 
respect to the higher modes interactions, it appears to be necessary to solve 
the first few modes simultaneously and to model the interactions from the rest 
of the modes. In that case, the difference between the first few modes in the 
role of forming a complete flow should become clear and this should answer the 
earlier question of the number of modes essential for the flow. The aniso- 
tropy in the interactions involving the higher modes will also decrease since 
the higher modes are expected to be more isotropic than the lower modes. 
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Table 1 Comparison of Reynolds stresses based on two sets of 
parameters: S = 0.01 and v T = 1 8 (Case I); S = 0.03 and 

v T = 22 (Case II). 
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Fig, 1 Turbulent kinetic energy versus t for various values of skewness 
parameter, S 
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Fig. 3 Profiles of 
numbers at t = 0.2 
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Fig. 5 Profiles of large eddy self- interact ions and modeled higher modes 
interactions versus y/d for = 5.0 and kg = 0.05 (1/cm) at t = 0.2 
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Fig. 9 Contribution from the first mode to (a) u , (b) v , (c) w 

— 2 

(d) uv, divided by U , versus y/d at t = 0.2: •• baufer [22] 




Fig. 10 Structural quantities calculated From the first mode versus y/d 
at t = 0.2: •• Laufe r [22]; (a) uv/q 2 , (b) uv//u^Vv'^, (c) 9 in equa- 
tion (21), and (d) v 2 /u 2 
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